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PREMCE  AND  SUHMAIg 

In  a  recent  article  R.  Wengert  suggested  a  technique  for  machine 
evaluation  of  the  partial  derivatives  of  a  function  given  in  analyti¬ 
cal  form.  In  solving  nonlinear  boundary-value  problems  using 
quasilinearization  many  partial  derivatives  must  be  formed  analytically 
and  then  evaluated  numerically.  Wengert' s  method  appears  very 
attractive  from  the  programming  viewpoint  and  permits  the  treatment 
of  large  systems  of  differential  equations  which  might  not  otherwise 
be  undertaken. 

In  this  Memorandum  we  show  how  to  apply  the  technique  to  some 
problems  of  orbit  determination,  though  our  ultimate  goal  is  to 
handle  much  more  complex  problems,  such  as  arise  in  adaptive  control. 
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I.  INTRODUCTION 


The  method  of  quasilinearization  for  the  numerical  solution  of 

nonlinear  multi-point  boundary^alue  problons  is  discussed  in  Refs. 

1-3.  Applications  to  system  identification,  cardiology  and  design 

and  control  are  given  in  Refs.  4-6.  The  determination  of  orbits  is 

discussed  in  Ref.  1.  The  method  involves  the  calculation  of  many 

partial  derivatives  which  can  be  quite  onerous  the  source  of  er- 

n\ 

rors.  In  a  recent  note'  R.  Wengert  suggested  a  procedure  for  the 
calculation  of  partial  derivatives  which  relieves  the  analyst  and  the 
programmer  of  much  routine  differentiation.  The  purpose  of  this 
Memorandum  is  to  show  the  usefulness  of  the  Wengert  scheme  in  carrying 
out  the  quasilinearization  procedure.  Our  discussion  is  built  around 
the  orbit  determination  problem  discussed  in  Ref.  1. 
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II.  ORBIT  DETEBMINATION 

Consider  a  heavenly  body  H  moving  about  the  sun  in  a  planar  or¬ 
bit,  the  x-y  plane.  The  equations  of  motion  are 


II 

X 


-X 


/  2^  2. 
(x  +  y  ) 


1.5 


f(x,y). 


(1) 


It 

y  = 


JiL. 


/  2_^  2, 

(x  +  y  ) 


g(x,y) 


(2) 


At  various  times  t^,  i=l,2,***,N,  the  angle  between  the  vector  from  a 

fixed  observer  at  (1,0)  and  the  x-axis  is  measured,  tan  e(t.)  ^  b, , 

1  i 

which  results  in  the  multi-point  boundary  conditions 

b^  =  tan  0(t^)  =  y(tp  /(x(tp-l),  i=l,2,***,N.  (3) 

We  wish  to  determine  a  set  of  conditions  on  x,  x,  y  and  y  at  an 
"initial  instant"  t=tj^  so  that  the  solution  of  Equations  (1)  and  (2), 
subject  to  these  conditions,  agrees  as  well  as  possible  with  the  ob¬ 
servations  in  the  sense  of  the  method  of  least  squares;  i.e.,  we  wish 
to  minlinlze  the  sum  S,  where 


2 

S  -  I  {y(tj)  -  [x(l:p  -l]  }  . 
i-1 


(4) 
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III.  QUASILINKARlZmON 

(1.2,4) 

Our  computational  formalism  proceeds  as  follows; 

We  first  select  an  initial  approximation  to  the  initial  conditions  and 
integrate  the  Equations  (1)  and  (2)  numerically  to  produce  the  func¬ 
tion  X  (t)  and  y  (t)  on  the  interval  t  ^  t  i  t^.  Then  we  proceed 

O  O  in 

recursively.  The  (k+1)®^  approximation  is  produced  from  the  k  by 
solving  the  linear  multi-point  boundary-value  problem 


+  (^k+r  V  yk>  ■ 


(5) 


(6) 


N 


Min 


k+1  ^*^i^  ’’  ^i  [\+l^*^i^  ’  } 


i=l 


(7) 


♦  • 

where,  the  minimization  is  over  x^(tj^),  x^^Ctj^), 

Particular  and  homogeneous  solutions  of  Equations  (5)  and  (6)  are 
readily  produced  numerically,  and  the  initial  conditions  are  chosen 
so  as  to  minimize  the  sum  in  Equation  (7).  See  Ref.  2. 
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IV.  WENGERI'S  METHOD 

Forming  the  partial  derivatives  called  for  in  Equations  (^5)  and 
(6)  is  not  especially  formidable.  If  other  sources  of  the  gravita¬ 
tional  field  were  present,  this  would  not  be  the  case.  Let  us  now 
see  how  we  can  have  the  computing  machine  evaluate  the  needed  partial 
derivatives . 

t-j\ 

Wengert's  method'  ^  is  based  on  the  chain  rule, 


du  _  3u  ^1  Bu  ^2 

dt  dt  Bx^  dt 


_  dx 
,  3u _ n 

cK  dt  ’ 
n 


(8) 


where  u  is  a  function  of  x,,  x„,  ••*,  x  ,  and  x  ,  x«,  •••,  x  are  func- 

iC  XX  jL  jL  xx 

tions  of  t.  By  putting  dx^^/dt  =  1  and  dx^/dt  =  0  for  i=2,*»»,N,  we 
find,  for  example,  du/BXj^  =  du/dt.  In  addition,  the  technique  uses 
the  elementary  formulas  of  differentiation.  The  reader  is  referred  to 
the  paper  by  Wengert  and  to  that  by  Wilkins.'  The  elementary  dif¬ 
ferentiation  subroutines  INTEXP,  ADD,  MCNST,  and  DIV  used  in  the  pre¬ 
sent  orbit  determination  program  have  been  described  in  Ref.  8.  These 
are  the  integer  exponentiation,  addition,  multiplication  by  a  constant, 
and  division  routines.  In  addition  to  these,  we  have  written  subroutine 
POW  for  non-integer  exponentiation, 


z  =  X 


z  =  c  X 


c-1  . 

X 


as  here  listed. 


SUBROUTINE  POW(X,C,Z) 

DIMENSION  X(2),  Z(2) 

Z(l)  =  X(1)**C  (z  = 

Z(2)  =  C*X(1)**(G-1.0)*X(2)  (z  =  c  x'^"^  x) 

RETURN 

END 
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Let  us  now  combine  Wengert*s  method  with  quasilinearization  and 
incorporate  the  resultant  scheme  into  the  FORTRAN  orbit  program.  But 
first,  let  us  write  the  linear  differential  equations  for  the  particu¬ 
lar  and  homogeneous  solutions  of  Equations  (5)  and  (6).  Let  p(t)  be 
the  particular  solution  corresponding  to  the  variable  x(t),  and  let 

q(t)  be  that  for  the  function  y(t).  Let  the  four  nomogeneous  solutions 

12  3  4 

for  x(t)  be  represented  by  h  (t),  h  (t),  h  (t),  and  h  (t) ,  Similarly 

12  3 

for  y(t),  we  have  the  homogeneous  solutions  w  (t),  w  (t),  w  (t),  and 
4 

w  (t).  Then  the  equations  for  the  particular  and  homogeneous  solutions 
are 

p  =  £(x^,  y^)  +  (P  -  y^)  +  (1  -  y^)  y^)  C) 

q  =  y^)  +  (P  -  V  V  +  -  yfc> 

h’’  =  y|j)  +''^£y<\-  y^)>  (11) 

»’■  =  y,^)  +  y,^).  (12) 

i  =  1,  2,  3,  4. 

Equations  (9),  (10),  (11),  (12)  show  that  the  functions  f,  f,, 
f  >  8>  8  >  S  must  be  evaluated  at  each  integration  step.  We  write 

y  ^  y  is; 

2  2  ^ 

subroutine  FUNl  to  evaluate  the  function  f(x,y)  =  -x/(x  +  y  )  ,  by 

calling  the  appropriate  elementary  differentiation  subroutines.  When¬ 
ever  this  subroutine  is  executed,  the  derivative  of  f  is  simultaneously 
computed.  The  following  is  a  listing  of  FUNl.  Beside  each  FORTRAN 
statement,  we  have  added  the  equivalent  mathematical  expressions. 
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SUBROUTINE  FUN1(X,Y,MS) 


DIMENSION  X(2) ,  Y(2) ,  ANS(2) ,  Zl(2) ,  Z2(2) ,  Z3(2) 
CALL  INTEXP(X,2,Z1)  Zj^  = 

CALL  INTEXP(Y,2,Z2)  ’ 

CALL  ADD(Z1,Z2,Z3) 

CALL  P0W(Z3,1.5,Z1) 

CALL  MCNST(-1.0,X,Z2) 

CALL  DIV(Z2,Z1,ANS) 

RETURN 


Z2  =  y 

^3  =^2 


h  =  (=^3^ 

22  =  -X 

f  =  Zj/'l 


1.5 


END 


z^  =  2xx 

22  =  2yy 

23  »  2i  +  22 

z^  =  Z3 

i2-i 

f  «  (22-fz.)/2j, 


Similarly,  we  write  subroutine  FUN2(X,Y,ANS)  to  evaluate  the  function 

g(x,y)  =  -  y/(x^+  y^)^*^. 

Next,  we  write  the  corresponding  partial  derivative  evaluation 
subroutines.  For  example,  the  following  subroutine  PDR],(X,Y,PD)  com¬ 
putes  f^  and  f^  and  stores  these  values  in  PD(1)  and  PD(2),  respec¬ 
tively  . 


SUBROUTINE  PDR1(X,Y,PD) 
DIMENSION  X(2) ,  PD(2) ,  A(2) 
X(2)  =1.0 
Y(2)  =0.0 
CALL  FUN1(X,Y,A) 

PD(1)  »  A(2) 

X(2)  =0.0 
Y(2)  -  1.0 
CALL  F0N1(X,Y,A) 

?n(2)  -  A(2) 

RETURN 

ENl^ 


X  »  1 
y  -  0 

f^  »  f,  since  X  «  1,  y  «  0 


X  -  0 
y  -  1 

fy  ■  f,  since  x  «  0,  y  ■  1 
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We 


write  also  the  partial  derivative  routine  for  and  g^,  subroutine 


PDR2. 

The  major  steps  involved  in  producing  the  right-hand  sides  of  the 
differential  Equations  (9)  to  (12;  are 


1)  CALL  EUNl  to  produce  yj^) 

2)  CALL  PDRl  to  produce  f  ,  f 

X  y 

3)  Evaluate  p,  h^,  h^ 

4)  CALL  FDN2  to  produce  g(Xj^»  y^^) 

5)  CALL  PDR2  to  produce  g^, 

6)  Evaluate  q,  w^,  w^,  w^,  w^. 
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V.  NDMERICAL  RESULTS 

We  repeated  one  of  the  numerical  experiments  mentioned  in  our  ear- 
lier  paper'’  * .  The  1962  experiment  was  carried  out  on  the  IBM  7090  by 
means  of  a  FORTRAN  II  source  program.  This  time,  using  Wengert's  me¬ 
thod,  we  ran  the  problem  on  the  IBM  7044  with  a  FORTRAN  IV  program. 

The  initial  approximation  to  the  orbit  is  a  fixed  point  over  all  time, 
this  point  coinciding  with  the  position  of  the  earth.  Seven  iterations 
were  carried  out,  and  the  results  in  both  runs  converged  to  the  correct 
solution  as  shown  in  Table  1. 


Table  1 


1962 

Experiment 

1964 

Experiment 

True 

Values 

x(2.5) 

1.19358 

1.19358 

1.19361 

x(2.5) 

-.664268 

-.664267 

-.664263 

y(2.5) 

1.06063 

1.06063 

1.06070 

y(2.5) 

.247466 

.247466 

.247499 

The  computer  execution  times  involved  are  1  minute  20  seconds  in  1962, 
and  2  minutes  35  seconds  in  1964.  It  is  difficult  to  estimate  the  in¬ 
crease  in  computing  time  with  the  use  of  Wengert's  method,  since  the 
computing  machines  and  the  source  languages  differed  in  these  two  rms. 
In  the  intermediate  iterations  1  through  5,  the  predicted  values  of  x, 
X,  y  and  y  at  time  2.5  varied  from  6  figures  of  agreement  between  runs. 
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to  only  2  figures  of  agreement.  The  final  results,  however,  did  agroe 
to  at  least  5  significant  figures. 

Die  use  of  Wengert’s  method  in  quasilinearization  calculations 
appears  promising,  and  additional  e3sperlments  will  be  carried  out. 
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